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We resolve a longstanding puzzle for the static and dynamic critical behavior of Gadolinium 
by a combined theoretical and experimental investigation. It is shown that the spin dynamics of 
a three dimensional ferromagnet with hep lattice structure and a spin-spin interaction given by 
both exchange and dipole-dipole interaction belongs to a new dynamic universality class, model J*. 
Comparing results from mode coupling theory with results from three different hyperfine interaction 
probes we find quantitative agreement. The crossover scenario for the wavevector dependence of 
the hyperfine relaxation rate is determined by a subtle interplay between three length scales: the 
correlation length, the dipolar and the uniaxial wave vector. 

PACS numbers: 75.40.Gb, 75.40.Cx, 76.75.+i, 76.80.+y 
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The spin dynamics of simple ferromagnets in the vicin- 
ity of their Curie point T c are archetypical examples 
of dynamic critical phenomena near second-order phase 
transitions. Much experimental and theoretical effort 
has been put into identifying the dynamic universality 
classes and assigning them to magnetic substances. Nev- 
ertheless, experimental observations on gadolinium [1-7] 
remained a puzzle up to now. Because of its large lo- 
calized magnetic moment, and the fact that it is an S- 
state ion, Gd should have a very small magnetocrystalline 
anisotropy and therefore be much better a model system 
for an isotropic Heisenberg magnet than either Fe, Ni or 
EuO. As a consequence it should belong to the model J 
dynamic universality class in the classification scheme of 
Ref. [8]. The measured static and especially dynamic 
critical exponents are, however, not at all compatible 
with model J. The objective of this paper is to resolve 
this longstanding seemingly contradictory situation by a 
combined theoretical and experimental study. 

Early experimental observations [9] clearly demon- 
strate that Gd has an easy axis which coincides with the 
hexagonal axis of its hep lattice. The origin of such an 
easy axis cannot be understood from the magnetocrys- 
talline anisotropy. But, based on a mean field theory [10] 
it has been argued that a combined effect of the lattice 
structure and dipolar interaction favors the c-axis as the 
easy direction. This view is supported by measurements 
[7] of the c-axis and basal-plane susceptibility on a single 
crystal of Gd. It is found that the basal-plane susceptibil- 
ity crosses over from a singular behavior to a constant at a 
characteristic temperature scale which can be accounted 
for by dipolar effects. The analysis of the static criti- 
cal behavior [7] is, however, complicated by the fact that 
all experiments are done in the non-asymptotic regime 
where superposed crossover lead to complex temperature 
dependences. This may not be easily interpreted in terms 
of one or the other universality class. This is even more 
so, as the static critical exponents for the various univer- 
sality classes are of comparable magnitude. 



A surprising and yet unexplained observation was 
made by a measurement of the critical dynamics us- 
ing hyperfine methods [1-3]. The critical exponent w 
for the autocorrelation time r c , which should scale as 
r c oc (T — T c )~ w in the asymptotic regime, was found to 
be w ~ 0.5. The observed value is not consistent with ei- 
ther Heisenberg or Ising models, but considerably lower. 

The purpose of this letter is twofold. First we give a 
theoretical description for a spin system with both ex- 
change and dipolar interaction on a hep lattice using 
mode coupling theory. Next we calculate the relaxation 
rates observed in various hyperfine interaction measure- 
ments, where we account for the details of the coupling 
tensor in each of these methods. We also report on mea- 
surements of the muon spin relaxation (/iSR) rate in high 
purity single crystal samples of Gd. A comparison of the 
theoretical predictions with these and earlier [4,6] /xSR 
measurements as well as perturbed angular correlation 
(PAC) and Moessbauer data [1-3] gives a coherent pic- 
ture of the dynamic and static critical behavior of Gd 
and resolves the puzzling situation described above. 

We consider a system with N identical spins fixed on 
the sites of a three dimensional lattice. Taking into ac- 
count magnetocrystalline anisotropy as well as dipolar 
interaction it is described by the Hamiltonian 
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The magnitude of the magnetocrystalline anisotropy of 
the system is characterized by A = j" /J ± . The dipolar 
interaction is characterized by the tensor 
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with Xy = Xj — Xj, <7l is the Lande factor, and hb 
the Bohr magneton. Dipolar lattice sums, = 

'e iq ' x * , can be evaluated by using Ewald's 
method. For infinite three-dimensional cubic lattices the 
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results may be found in Ref. [11,12]. For Bravais lat- 
tices with a hexagonal-closed packed (hep) structure the 
dipolar tensor to leading order in q becomes [10] 
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where v a is the volume of the primitive unit cell with 
lattice constant a, and the parameters are /3f = 4.12 
and /3| = 4.32. Upon expanding the Fourier trans- 
form of the exchange interaction = ^ J$e iq ' Xi w 
Jq — Jq 2 a 2 + 0(q 4 ), and keeping only those terms which 
are relevant in the spirit of the renormalization group 
theory, one finds 
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There are two sources of uniaxial anisotropy in the 
Hamiltonian, Eq. 4, magnetocrystalline anisotropy, Aq = 
Jq/J, and dipolar interaction, m a — (gLl^B) 2 P2 /2Jiv 
In addition, the dipolar interaction introduces an 
anisotropy of the spin-fluctuations with respect to the 
wave vector q which is reflected by the term propor- 
tional to q a q(j/q 2 . The magnitude g of this anisotropy is 
given by g = ^{gh^B) 2 /2Jv a . We define a dimension- 
less quantity m = {gL^B) 2 {fi\ — Pir)/2Jv a proportional 
to the ratio of the anisotropy energy and the exchange 
energy. Putting in values for Gd the ratio of the dipo- 
lar contribution to the term q a q(}/q 2 and to the uniaxial 
anisotropy is \Jgjra = 7.8738. In the following we will 
show that all available data for Gd can be explained by 
assuming that the uniaxial anisotropy is solely due to the 
dipolar interaction. Therefore, we will assume A = in 
the following discussion. 

Now we turn to an analysis of the static critical behav- 
ior. In Ornstein-Zernikc approximation the static suscep- 
tibility reads 
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where r z = r = £~ 2 , r x , y = S, 2 + Qa an d we have 
measured all length scales in units of the lattice con- 
stant a. The analysis of the critical behavior resulting 
from the Hamiltonian, Eq. 4, is complicated by the fact 
that besides the correlation length £ = £ (T/T C — 
there are two anisotropy length scales q^ 1 — a/y/m 
and q^ 1 = a/y/g, both resulting from the dipolar in- 
teraction. The eigenvalues of the inverse susceptibil- 
ity matrix are given by Ai(q) = q 2 + + q\ and 
A 2 , 3 (q) = q 2 + + [q 2 D +q 2 A ±W] /2 where W = 
{(q 2 D + q\) 2 - 4:q 2 D q 2 A q 2 /q 2 } 1/2 ; the eigenvectors ej(q) are 
given in a forthcoming publication [13]. It is interesting 
to note that due to the combined effect of the dipolar 
interaction and the uniaxial anisotropy of the lattice the 



eigenvalues of the susceptibility matrix remain finite in 
the limit q — > and upon approaching the critical tem- 
perature. Only if the angle •& between the easy axis (z- 
axis) of magnetization and the wave vector is $ = 90° the 
third eigenvalue becomes critical. Due to the anisotropy 
terms in the Hamiltonian there are two crossover as a 
function of the reduced temperature and wave vector. 
To a good approximation these can be incorporated in 
an effective exponent of the correlation length £ [14]. 

Mode coupling theory is a theoretical method which 
has been shown to give highly accurate results for the 
critical dynamics of cubic ferromagnets [12]. Here we 
generalize this method to non-cubic magnets. Starting 
from the equations of motion for the components Sq of 
the spin S q in the eigenvector basis, = J2i ^q e ai{<l), 
one can derive the following set of coupled integral equa- 
tions [13] for the half-sided Fourier transform, ^(q, w) = 
«Xa(q)/[w + «r Q (q, oj)], of the Kubo relaxation function 
$a/3(q,i), 



r a (q,w) = 
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where j*k,w = J(d 3 fc/(27r) 3 ) J(du}/2n). The vertex func- 
tions if^ 7 (k, q) for the decay of the mode a into the 
modes [3 and 7 are given by [13] 

i^(k, q) = (k, q)[&(k, q) [A^k) - A^q - k)] , 



K^(k, q) = Tf(k, q) T^(k, q), ? 7 , 



(7) 



with f/^(k,q) = ^ y . fe e 4jfe e cti (k)e / 3 i (q)e 7 fc(q - k), 
where e^k is the Levi-Cevita symbol, and T"(k, q) = 
A 4 (k)L^(k,q), T£(k,q) = [A i (k)-A,(q-k)]f/4(k,q) 
for i < j and T^(k, q) = for i > j. In the limits 
qo — > or qA — > these equations reduce either to the 
uniaxial or to the isotropic dipolar ferromagnet [12,15] 

We have solved the above mode coupling equations in 
the Lorentzian approximation for the lineshape [13]. It 
is found that the mode coupling equations obey a gener- 
alized dynamic scaling law where the linewidth depend 
on three scaling variables, x\ = l/g£, x^ — qo/q and 
x 3 = qA/q, and the angle $ between the direction of the 
wave vector and the z-&xis. The explicit functional form 
is too complicated to be presented in this letter and we 
refer the reader to a forthcoming publication [13]. But, 
note that it is implicitly contained in the damping rates 
for the hyperfine interaction probes discussed below. 

Now we compare the above theoretical results with ex- 
perimental data obtained from hyperfine interaction ex- 
periments. We have performed zero field /iSR experi- 
ments on Gd [16]. Spin polarized muons were implanted 
in order to measure the distribution and the dynamics 
of the internal magnetic field at the muon site via the 
temporal loss of the initial spin polarization. The exper- 
iments were performed at the ^SR facility of the Paul 
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Scherrer Institute using a high momentum muon beam. 
The sample was a spherical Gd single crystal with di- 
ameter 2.5 cm. The temperature could be stabilized to 
at least ±0.05K . We could describe the temporal loss 
of the initial muon spin polarization by an exponential 
decay function P(t) = cxp(— X z t). Taking into account 
anisotropic dipolar fields as well as the isotropic Fermi 
contact field to the local field at the muon site the muon 
damping rate can be written as [17,18,13] 



A* 



W q 



E 



[GfG% + GfG% 



$^(q,0), (8) 



where the hatted variables indicate that the correspond- 
ing quantities have to be evaluated in the muon ref- 
erence frame, i.e. the z-axis coincides with the initial 
polarization of the muon beam. Here we have defined 
T> = 7^(/j,o/47t) 2 ((7l^.b) 2 . The coupling of the muon spin 
and the spins of the magnet is described in terms of the 
coupling matrix G™' 3 , which reflects the particular sym- 
metry of the lattice sites occupied by the muons. Since 
the most dominant contribution to the damping rate 
comes from wave vectors close to the Brillouin zone cen- 
ter, the form of the coupling tensor G 3 ^ at small values of 
q will be important. In the limit q — > [18,13] one finds 



G. 



af3 



-47T [q a qp/q 2 - p a ] , with p a = d a + n^H^/Air. 
With the Fermi contact field B FC = -6.98 kG at T = K 
[19], one gets n^H^/An = —0.278 [18], and consequently 
for octahedral sites p x = p v = 0.0705 and p z = 0.0250 
[18], and for tctrahcdral sites p x = p y = 0.0338 and 
p z = 0.0984 [13]. The muon relaxation rate A z depends 
on the material parameters q^Co and <?_d£o- Since we 
assume that both anisotropies result from the dipolar 
interaction the ratio qo/qA = 7.8738 is known and the 
number of material parameters is reduced to one. In com- 
paring our theory with fiSK experiments at a polarization 
a = 90° we get the best fit to the data with qo£,o = 0.13 
(see Fig. 1). This results in the following values for 
the uniaxial and dipolar wave vector, q& = 0.0165/£oj 
qo = 0.13/£o- The corresponding crossover temperatures 
are given by T A =T C + 0.43 K and T D = T c + 16.54 K 
suggesting the following crossover scenario. For T 3> To 
we expect critical behavior dominated by the (isotropic) 
Heisenberg fixed point. The relaxation rate shows power 
law behavior A oc t~ w , with an exponent w w 1. For 
temperatures in the interval To > T > Ta dipolar in- 
teraction becomes important. But, from the analysis of 
the uniaxial crossover [13] it turns out that the uniax- 
ial crossovers in dynamics sets in at wave vectors much 
larger than expected from an analysis of the static quan- 
tities. Therefore, even for T > Ta we expect to ob- 
serve effects from dipolar interaction as well as uniaxial 
anisotropy. Finally, for T < Ta the critical dynamics is 
determined by the uniaxial dipolar fixed point. Then the 
static susceptibilities do no longer diverge for q — > and 
T — > T c except when the wave vector is perpendicular to 



the easy axis of magnetization. Since the relaxation rate 
A z is given by an integral over the whole Brillouin zone, 
the relative weight of the critical axis along which the 
susceptibility diverges becomes vanishingly small. As a 
consequence the relaxation rate A z no longer diverges for 
T — > T c , i.e. w — > (compare Fig. 1). 

Fig. 1 shows a comparison between the theoretical and 
experimental results for an initial polarization inclined 
by an angle a — 90° with respect to the easy axis of 
magnetization. The solid line and the dashed line are 
the theoretical result for the muon relaxation rate if the 
muons penetrating the sample are located at tetrahcdral 
and octahedral interstitial sites, respectively. The com- 
parison between theory and experiment favors tetrahc- 
dral sites. This is confirmed by /iSR experiments with 
the initial polarization along the easy axis of magnetiza- 
tion. The ratio A z (90°)/A z (0°) for T -> T c becomes 1.2 
and 0.7 for octahedral and tctrahcdral sites, respectively 
[16]. The experiment is closer to the latter value strongly 
suggesting that muons occupy octahedral sites within the 
Gd lattice. 
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FIG. 1. Experimental and theoretical results of the relax- 
ation rate A for tetrahedral and octahedral muon sites with 
a = 90°. Data taken from Ref. [6,4] and measured at the /iSR 
facility at the PSI in 1992 and 1994 (see inset) [16] 

The coupling tensor in PAC and Moessbauer measure- 
ments reduces to a Fermi contact coupling. As a conse- 
quence the observed relaxation rate is a sum over the 
eigenmodes r c = (k B T/3v a ) Y, a J Xa(q)/ r Q (q)- Im- 
portant information about the behavior of the auto- 
correlation time can be gained from a scaling analysis. 
An effective dynamical exponent z c g may defined by 
t c oc (T — T c )~ Waii with to e ff = v s{z c s — 1), where we 
have neglected the Fisher exponent rj. If dipolar inter- 
action and uniaxial anisotropy were absent, one would 
expect to w 1.0. The dipolar interaction is known to be 
a relevant perturbation at the Heisenberg fixed point. It 
leads to asymptotic static critical exponents which are 
only slightly different from the corresponding Heisenberg 
values. But, since dipolar interaction implies a non- 
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conserved order parameter the asymptotic dynamic ex- 
ponent becomes zd ~ 2 resulting in a crossover from 
Wi w 1.0 to wd ~ 0.7. Uniaxial interaction is also known 
to be a relevant perturbation with respect to the Heisen- 
bcrg fixed point. Again, the static critical exponents are 
not changed very much, e.g. one finds the Ising (I) value 
v\ = 0.63, but the dynamic exponent becomes z\ w 4 if 
the order parameter is conserved (z\ w 2 otherwise). The 
corresponding exponent for the hyperfine relaxation rate 
would turn out to be w\ 1.89 and Wi ss 0.63 for con- 
served and non-conserved order parameter, respectively. 
According to these scaling arguments it is hard to think 
of any dynamic universality class which could lead to an 
effective exponent w c s smaller than about 0.6. Actually, 
however Moessbauer studies and PAC measurements on 
Gd show distinctly anomalous low values w w 0.5, which 
cannot be explained by either of the above scenarios. 
This experimental puzzle can be resolved if one consid- 
ers the combined effect of dipolar interaction and uniaxial 
anisotropy. As we have seen in the above analysis of the 
static critical behavior of uniaxial dipolar ferromagnets, 
all eigenvalues of the susceptibility matrix remain finite 
upon approaching the critical temperature except when 
the wave vector of the spin fluctuations is perpendicular 
to the easy axis. Since this is only a region of measure 
zero in the Brillouin zone one actually expects that the 
relaxation rate does no longer diverge upon approaching 
T c , i.e., w — > 0. 



wave vector as for our comparison with ^iSR experiments. 
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FIG. 2. Experimental and theoretical results of the auto- 
correlation time t c for PAC-experiments. Data taken from 

[!]■ 
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FIG. 3. Experimental and theoretical results of the auto- 
correlation time t c for Moessbauer spectroscopy experiments. 
Data taken from [2,3]. 

In summary, we have outlined a mode coupling theory 
for uniaxial dipolar ferromagnets, where the uniaxiality 
solely results from the dipolar interaction. We have also 
reported measurements on a high purity single crystal 
sample of gadolinium metal in the paramagnetic regime. 
From the quantitative agreement between this theory and 
our fiSR and previous PAC and Moessbauer measure- 
ments the following conclusions can be drawn: (i) The 
universality class of Gd is the uniaxial dipolar ferromag- 
net, where both the isotropic dipolar and the uniaxial 
contribution to the spin Hamiltonian are due to a com- 
bined effect of dipolar interaction and non-cubic lattice 
structure. This corresponds to a new anisotropic model 
J*, where the Hamiltonian is given by Eq. 4. (ii) The 
dominant factor for the uniaxial anisotropy in Gd is the 
dipolar interaction, (hi) Muons in Gd are located at oc- 
trahedral interstitial sites close to T c . 
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Let us now compare the results of our mode cou- 
pling theory with hyperfine experiments on Gd men- 
tioned above [1-3] . The autocorrelation time r c is shown 
in Figs. 2 and 3 for PAC experiments and Moessbauer 
spectroscopy, respectively. Both sets of data are in ex- 
cellent agreement with the results from mode coupling 
theory for T - T c < 10K. Note that besides the overall 
frequency scale there is no fit-parameter, since we have 
used the same set of values for the dipolar and uniaxial 
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